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ABSTRACT 

We study the structure of passively heated disks around T Tauri and Herbig Ae stars, and present a vectorized Monte 
Carlo dust radiative transfer model of protoplanetary disks. The vectorization provides a speed up factor of ~ 100 when 
compared to a scalar version of the code. Disks are composed of either fluffy carbon and silicate grains of various sizes 
or dust of the diffuse ISM. The IR emission and the midplane temperature derived by the MC method differ from 
models where the radiative transfer is solved in slab geometry of small ring segments. In the MC treatment, dusty halos 
above the disks are considered. Halos lead to an enhanced IR emission and warmer midplane temperature than do pure 
disks. Under the assumption of hydrostatic equilibrium we find that the disk in the inner rim puffs up, followed by a 
shadowed region. The shadow reduces the temperature of the midplane and decreases the height of the extinction layer 
of the disk. It can be seen as a gap in the disk unless the surface is again exposed to direct stellar radiation. There 
the disk puffs up a second time, a third time and so forth. Therefore several gaps and ring-like structures are present 
in the disk surface and appear in emission images. They result from shadows in the disks and are present without the 
need to postulate the existence of any companion or planet. As compared to Herbig Ae stars, such gaps and ring-like 
structures are more pronounced in regions of terrestrial planets around T Tauri stars. 
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1. Introduction 

Low-mass stars are surrounded by gaseous dust disks; 
prominent examples are T Tauri and Herbig Ae/Be stars 
(Waters & Waelkens, 1998). The existence of disks is de- 
rived from observations of the submillimeter continuum 
(Beckwith et al., 1990), the IR excess over the photospheric 
emission (Kenyon & Hartmann, 1995), and direct imaging 
(McCaughrean & O'dell, 1996). Protoplanetary disks with 
weak mid-IR and strong far-IR emission are called transi- 
tional disks (Najita et al. 2007; Furlan et al. 2009; Cieza et 
al. 2008). They are of interest because the missing mid-IR 
emission could be caused by grain growth or planet for- 
mation (Testi et al., 2003). Their optical extinction is flat- 
ter than the standard ISM extinction curve (Fitzpatrick & 
Massa, 2007) and shows time variability (Schisano et al. 
2009). This indicates that large grains or clumpy material 
orbits and partially obscures the star. 

The evolution of the disk is followed theoretically by 
studying the time dependence of the surface density of gas 
and dust. Hydro-dynamical models of viscous protoplane- 
tary disks including planets have been developed (Lubow et 
al. 1999; Ozernoy et al. 2000; Klahr & Bodenheimer 2003; 
Edgar & Quillen 2008; Alexander & Armitage 2009). The 
models predict that the disk is cleared by an orbiting planet 
and that the surface density shows structures with holes 
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and gaps, spiral arms, and rings. Such disk structures are 
inferred from fitting the spectral energy distribution (SED) 
by applying models considering pure disks (Calvet et al 
2002; Espaillat et al. 2010; Mulders et al. 2010) or disks 
with halos (Schegerer et al. 2008; Verhoeff et al. 2011). 

Recently, several disks with rings have been observed 
(van Boekel et al. 2005; Fukagawa et al. 2006; Ratzka et 
al. 2007; Brown et al. 2007, 2008; Mayama et al. 2010; 
Thalmann et al. 2010; Olofsson et al. 2011). In the gaps, 
between the rings, companion candidates are detected 
(Huelamo et al. 2011; Kalas et al. 2008; Lagrange et al. 
2010). Hydro- models account for the observed radial distri- 
bution of planets and lifetimes of the disks, which are a few 
Myr (Strom et al. 1989; Bertout et al. 2007). The strength 
of the disk-planet interactions are influenced by the gas 
pressure so that the structure of protoplanetary disks plays 
a key role in understanding the planet formation process. 

There is a rich literature on radiative transfer models 
of protoplanetary disks (e.g. Wolf et al. 2003; Whitney et 
al. 2003; Pinte et al. 2006; Robitaille et al. 2006; Pintc 
et al. 2008). The structure of the disks is discussed by 
Nomura (2002), Walker et al. (2004), and Dullemond & 
Dominik (2004a), and self-shadowing ripples in the disks 
have already been mentioned by D'Alessio et al. (1999) and 
Dullemond (2000); for a recent review, see Dullemond & 
Monnier (2010). In this paper the thermal structure of disks 
is discussed by solving the hydrostatic and radiative balance 
equations. We are interested in deriving the structure and 
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IR appearance of a starlight -heated, so-called passive disk. 
Dust models used in this paper are described in Sect. [2j 
and a vectorized three-dimensional Monte Carlo (MC) ra- 
diative transfer code in Sect. [3] We compare results of the 
MC model against a method that samples the disk in small 
ring segments at radius r from the star and in which the 
radiative transfer of the segments is solved in a slab ge- 
ometry (Sect. |4~3")) . The influence of different dust opacities 
on the SED of the disks is discussed in Sect. 14.41 We ap- 
ply the models to T Tauri and Herbig Ae stars. The effect 
of a dust halo on the emission and temperature structure 
is presented (Sect. l4"3j) . The vertical temperature structure 
introduces variations of the height of the disk surface which 
causes shadows. They emerge as gaps and ring-like struc- 
tures at distances from the star, which are thought to be 
the birthplaces of terrestrial planets (Sect. [SJ. Strikingly 
the structures are caused only by the coupling of the hy- 
drostatic and radiative transfer equations without the need 
to postulate a planet. 

2. Dust models 

In protoplanetary disks the properties of dust are thought 
to evolve, and the mineralogy, composition, porosity and 
size distribution of the grains differs from that of the diffuse 
ISM. Thanks to an enhanced particle collision rate, fluffy 
grains are produced that grow in size over time (Natta et 
al. 2004; Acke et al. 2004; Natta et al. 2007; Lommen et al. 
2007; Ricci et al. 2010). Disk particles settle towards the 
midplane under gravity, a process counteracted by grain 
fragmentation and disk turbulence (Mizuno et al. 1988; 
Weidenschilling et al. 1993; Sterzik et al. 1994; Dullemond 
& Dominik 2004b, 2005). It is found that grains of radius 
>10jUm settle towards the midplane, and depending on tur- 
bulent motions, will be blown up to the surface, where they 
remain for long time periods, and well coupled with the gas 
(Alexander & Armitage 2007; Charnoz et al. 2011). 

Dust gets annealed in high-temperature regions of the 
disk, and crystalline structures are built. They are detected 
and the observed profiles of the 10//m silicate band differ 
from that of the interstellar medium (Malfait et al. 1998; 
Bouwman et al. 2008; Schegerer et al. 2006; Kessler-Silacci 
et al. 2006; Furlan et al. 2011; Oliveira et al. 2010; Watson 
et al. 2009). About 4% to 7.5% of the dust mass of the Solar 
System is in crystalline silicates whereas crystallization is 
~ 2% in the diffuse ISM (Kemper et al., 2005). Ice coatings 
are built in frosty regions (Siebenmorgen & Gredel 1997; 
Pontoppidan et al. 2005; Visser et al. 2009). The location 
where water ice can exist is important for terrestrial planet 
formation (Ida & Lin 2008; Min et al. 2011). All of these 
processes are combined by mixing and transport mecha- 
nisms and are valid for the neutral part of the disks where 
grain charging can be ignored. Ionization becomes impor- 
tant in the disk surface; therefore, the true situation of the 
mineralogy in protoplanetary disks is a complicated matter 
with a rich chemical network and evolution mechanisms in 
place (Gail, 2004). 

The SED and the thermal structure of the disks depend 
on the applied dust model and for comparison we present 



results using dust opacities as derived for protoplanetary 
disks and the diffuse ISM. For the dust in the ISM, we use 
a power-law size distribution: n(a) oc a -3 ' 5 with particle 
radii between a_ < a < a + of bare spherical particles of 
silicates (a_ = 320A , a+ = 2600A) with optical constants 
by Draine (2003) and carbon (a_ = 160A, a+ = 1300A), 
with optical constants by Zubko et al. (1996), and both with 
bulk density of 2.5g/cm 3 . We use dust abundances (ppm) of 
31 [Si]/[H] and 200 [C]/[H], which are in reasonable agree- 
ment with cosmic abundance constraints (Asplund et al., 
2009). This gives a gas-to-dust mass ratio of 130. For pro- 
toplanetary dust we consider fluffy agglomerates of silicate 
and carbon as subparticles. Other parameters are as for the 
ISM dust except that grains grow to 100 times larger radii 
of a+ = 33/im. As relative volume fraction of the composite 
grain we use 34% silicates, 16% carbon, and 50% vacuum, 
which translates to a relative mass fraction of 68% silicates 
and 32% aC. Absorption and scattering cross-sections and 
the scattering asymmetry factor is computed with Mie the- 
ory for the ISM grains, and we apply the Bruggeman rule 
for the fluffy composites. This gives a total mass extinction 
cross section for the large grains in the optical (0.55/im) 
of K^ xt = 30600 of the ISM dust and 4000 [cm 2 /g-dust] 
of the fluffy composites. The wavelength dependence K ? 
of both models is similar to those displayed in Kriigel & 
Siebenmorgen (1994, Fig. 12). 



3. Vectorized MC radiative transfer 

The radiative transfer problem is solved by a MC technique 
by following the flight path of many random photons. Basic 
ideas are given by Witt (1977), Lucy (1999), and Bjorkman 
& Wood (2001), and the original version of our code was 
developed by Kriigel (2008). Our model space is a three- 
dimensional Cartesian grid (x, y, z) that is partitioned into 
cubes and, where a finer grid is required, further divided 
into subcubes (cells) . The star of luminosity L emits a to- 
tal of N = n ■ N v photon packages per second and in each 
of the N„ frequency bins n photon packages are emitted. 
Each package has a constant energy e = L/N. A package 
entering a cell on its flight path may be absorbed there or 
scattered. The probability of such an interaction is given 
by the extinction optical depth along the path within the 
cell, At, and occurs when At > —log(£) using unified ran- 
dom number £. The package is scattered if the dust albedo 
A > using random number otherwise, it is absorbed. 
When the package is scattered, it only changes direction de- 
termined in a probabilistic manner by the phase function. 
When it is absorbed, a new package of the same energy, 
but usually different frequency v' is emitted from the spot 
of absorption. The emission is isotropic. Each absorption 
event raises the energy of the cell by e, and accordingly its 
temperature. 

The computational speed of the code is increased con- 
siderably by calculating the flight paths of photons simulta- 
neously. This type of vectorization is realized in two flavors 
using either conventional computer processing units (CPU) 
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within FORTRAN 90 and OMP environment!]] or graphics 
processing units (GPU) with NIVIDA cards and CUDAB 
(Heymann, 2010). The speed-up scales almost linearly with 
the number of processing cores available. For parallelization 
particular attention needs to be paid to the random num- 
ber generator, for which we choose the Mersene Twister 
algorithm (Matsumoto & Nishimura, 1998). 

Verification of the vectorized MC code against bench- 
marks (Ivezic et al. 1997; Pascucci et al. 2004) was done in 
an accompanying paper (Heymann & Siebenmorgen, 2012). 
The code is also tested against the spheric symmetrical ray 
tracing code by Kriigel (2008). The SED is computed by 
counting the photon packets leaving the cloud towards the 
observer using a large beam. The observed intensity can 
also be computed by following the radiative transfer of the 
line of sight from the observer or a pixel of the detector 
plane, through the model cloud using the MC computed 
dust temperatures and scattering events (Heymann, 2010). 
The ray tracer is also vectorized and allows us to derive dust 
scattering and emission images with high signal to noise, 
which is not possible by photon counting procedures. 

The number of interactions of photon packets with the 
dust increases exponentially with the optical depth of the 
cell. In cells of extreme high optical depth, ry > 1000, the 
photon packages are trapped and interact with the dust 
over and over again before they have a chance to escape. 
In this way MC treatments are slowed down considerably. 
Unfortunately, this situation appears in cells close to the 
midplane of protoplanetary disks. A modified random walk 
(MRW) procedure for improving the computational effi- 
ciency of MC methods has been presented by Fleck & 
Canfield (1984). The MRW enlarges the mean free path 
length of packets by a diffusion approximation whenever 
necessary, and has been tested by Min et al. (2009) and 
Robitaille (2010). Let r be the distance of the interaction 
point to the nearest site wall of a cell. The MRW is con- 
sidered in a cell when the photon package has a distance 
r > l/pKu, where 7 > 5 is a user specified constant, p 
the dust density, and Kp the Rosseland mean extinction 
coefficient. The photon travel distance d is derived from 

d=-^^)\s (1) 

with diffusion constant D = 1/SpKp and a pre-computed 
sample of S given by 

00 

£ = 2£(-iy +1 (<p) 2 (2) 

1=1 

with unified random number £. The number of absorption 
events within a cell is used to compute the dust temperature 
and, within r of the MRW, is estimated by n a b s — dp Kp, 
where Kp is the Planck mean mass extinction coefficient 
of the dust. Ignorance of these acceleration methods has a 
tremendous effect on the run time requirements of a vector- 
ized MC treatment. In our scheme the trajectory of photon 

1 www.openmp.org 
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packages through the model space is vectorized in so-called 
threads. Trajectories that hit the cells of high optical depth 
have many more interactions than others, most likely the 
majority of trajectories. This results in a rather unbalanced 
workload over all threads, so that the advantage of vector- 
ization is lost. In such cases the run time increases by the 
number of CPU or GPU cores available, which in our case 
are 8 and 480, respectively. This factor of efficiency loss is 
on top of the speed up gain of 5 - 10 of the MRW procedure 
achieved in scalar MC treatments (Pinte et al., 2009). 

4. Protoplanetary disk models 

The time scale for disks to attain thermal equilibrium is set 
by the balance between heating and cooling. This time scale 
is much shorter than the evolution of the disk or the heating 
source. The number density of the gas, even in the upper 
layers of the disk, exceeds 10 5 cm~ 3 , so that gas and dust 
are thermally coupled. Matter will spiral in the disk towards 
the star and dissipate gravitational energy. For a thin disk 
that extends to the stellar surface, an accretion luminosity 
of about a quarter of the stellar luminosity is converted in 
this way into heat. For classical T Tau stars, observed ac- 
cretion rates are ~ 10~ 7 '"~ 9 M Q /yr (Muzerolle et al., 1998). 
At early epochs when the accretion is strongest, dissipation 
dominates, but later heating by stellar radiation becomes 
more important. As fiducial model of the T Tau star, we 
take a mass of IMq, a luminosity of 2Lq, and a photo- 
spheric temperature of 4000 K. For comparison we treat a 
disk heated by a Herbig Ae star with 2.5M Q , 50L© and 
10 4 K (van den Ancker et al., 1997). 

^From near-IR interferometry (Millan-Gabet et al., 
2006) of T Tauri and Herbig Ae stars, it has been found 
that the inner disk scales with stellar luminosity n n oc \[L. 
Such a dependency is expected assuming that r; n is set 
by evaporation of grains at temperatures of 1000 — 1500 K 
(Dullemond & Monnier, 2010). For Herbig Be stars, larger 
deviations of the simple relation are measured at L > 
10 3 L o . Spatially resolved spectroscopy, at milliarcsec reso- 
lution, allows detection of a hotter gaseous emission compo- 
nent closer to the star (Eisner et al. 2009; Najita et al. 2009; 
Pontoppidan et al., 2011). Evaporation temperatures of sil- 
icate and carbon particles are around 1500 K. The grains in 
the disk are assumed to be fluffy and to have been formed 
by coagulation of much smaller compact interstellar grains. 
The fluffy agglomerates are held together by weak van der 
Waals forces (binding energies <0.1eV), and they are ex- 
pected to fall apart into their refractory constituents when 
the temperature of the composite grain exceeds 1000 K. The 
constituent particles will then be hotter than the average 
porous grains because they are much smaller. We assume 
that the disk extends inwards up to the point where porous 
grains reach 1000 K, or equivalently small interstellar grains 
would be at about 1500 K. 

In star-forming regions dust is heated to about 30 K. At 
a large distance from the star grain heating by the ambi- 
ent radiation field of nearby low or massive stars becomes 
important. We do not make additional assumptions on the 
outer radiation field and consider only the primary central 
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Table 1. Parameters of the fiducial T Tauri and Herbig Ae disks. 



Parameter 


T Tauri 


Herbig Ae 


Stellar luminosity L* [LqJ 
Stellar mass M* [Mq] 
Photospheric temperature T» [K] 

Column density E(r) - T ^ AU) ( A u) 7 [g dust/cm 2 ] 

Vertical optical depth rx(lAU) 

Dust density in halo Phaio [g-dust/cm 3 ] 

Inner disk radius n n 

Outer disk radius r out [AU] 


2 
1 

4,000 
r < 1AU 
r > 1AU 
1C 

or 1. 
evap 

22.5 


50 
2.5 

10,000 
: 7 = 0.5 
: 7 = -l 
,000 

5 x 10" 18 
oration 
40 



heating source. Models are computed to an outer radius of 
r out = 22.5 AU for T Tau and 40 AU for Herbig Ae disks 
where the midplane temperature drops below 30 K. 

Initially we assume that the disk is isothermal in z, and 
the density is given by 

with scale height H 2 = kT^^r 3 /GM^m, surface density 
S(r), molecular mass m = 2.3m p and midplane tempera- 
ture, T m id, for which we use a power law as initial guess. 
The height of the disk is set to z$ — 4.5 H. In models 
with pure disks, the density phaio above z > z a is and 
constant when a halo is considered (Table [1} . The surface 
density is adjusted to reach a given optical depth in the 
vertical direction t± = ~E(r)Ky, with visual extinction Ky. 
We use t_l = tiau (r/AU) 7 with 7 = — 1 for r > lAU and 
7 = 0.5 otherwise (Min et al. 2011). For the porous grain 
model with a+ = 33^m, we take a vertical optical depth 
from the surface to the midplane at 1 AU of t± = 10,000. 
This choice holds computational time requirements of the 
MC code within reasonable limits. It translates to a surface 
density of S(1AU) = 5 g-dust/cm 2 or a gas surface density, 
which is half the value estimated for the minimum mass 
of the early solar nebulae (Hayashi, 1981). Higher surface 
densities can be obtained considering larger or ice-coated 
grains because Ky is reduced in such models (Kriigel & 
Siebenmorgen, 1994). 

4.1. Disks with slab geometry 

A solution of the radiative transfer of a hydrostatic and 
geometrically thin disk is developed by Kriigel (2008, 
Sect. 11.3). The disk is symmetric with respect to the mid- 
plane at z = 0. The density structure is given in cylindrical 
coordinates p(r,z). Light from the star falls on the disk 
under a grazing angle a gr so that the star is visible from 
everywhere on the disk surface (Armitage, 2007). Flat disks 
with a constant grazing angle of cv gr = 2° and flared disks, 
similar to Chiang & Goldreich (1997), are considered with 
3° < cv gr oc r 2 / 7 < 7°. The flat disk has a smaller grazing 
angle and intercepts less stellar radiation than the flared 
disk, and even more so at large distances from the star. 



The disk is divided into small ring segments of width 
Ar at radius r from the star. For each ring the radiative 
transfer is solved under all inclination angles 9, for incom- 
ing I~ and outgoing I + radiation, assuming a slab geome- 
try. In the vertical direction the opacities are so high, e.g. 
tl (1AU) = 10 4 , that each ring segment is split into a com- 
pletely opaque midlayer and a much thinner top layer of op- 
tical thickness T top . The midlayer is assumed to be isother- 
mal at temperature T m ;d and the top layer extends so far 
down that the temperature at its bottom approaches T m id- 
For computing the spectral energy distribution it is suffi- 
cient to choose Tt p ~ 20. The ray tracer is an iteration 
procedure and yields the temperature structure T(r, z) and 
the intensities of the upward / + and downward / _ directed 
radiation. By observing the disk at viewing angle # bs> the 
received flux is computed from the emission of the star, 
the intensity I + at the disk surface at 8 = 8 b s as well as 
by summing up contributions from all rings. The radiative 
transfer is solved for each ring segment at radius r sepa- 
rately. However, the propagation of the radiation in radial 
direction from one ring into the next is ignored. The pro- 
cedure is often called the 1+1D disk model. 

4.2. Disks with axial symmetry 

In disks that are in hydrostatic equilibrium in a vertical 
direction, the gravitational force is balanced by the pressure 
gradieno 

z GM* _ 1 dP 
tt 2 p dz 

with pressure P = pkT(z)jm. For each height z and ra- 
dius r 2 = x 2 + y 2 from the star, the temperature T(r, z) 
is computed using the MC code of Sect. [3l The star is set 
at origin. Azimuthal dependencies in the disk structure are 
not considered, so that the problem can be solved in axial 
symmetry. It allows us to solve the radiative transfer in one 
octant of a cube and reduces the computational time by a 
factor 8 when compared to a full three dimensional treat- 
ment. For each height z > and radial distance at x > 
and y > 0, we compute the azimuthal average of absorbed 

3 The protoplanetary disk masses are well below the stellar 
mass A/disk < 0.01M*. 
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photon packets. This average is used to derive the dust 
temperature T(r, z), and with Eq. [4] a new estimate of the 
density structure of the disk p(r,z). Obviously there is an 
iterative scheme: starting with an initial guess of the mid- 
plane temperature and a first set up of the disk structure 
p(r, z) using Eq. [31 we compute the temperatures T(r, z) 
after a MC run. Once the temperatures are inserted into 
Eq. 21 we arrive at a new density structure p(r,z), which 
is used to update T(r, z) with a new MC run. Typically 
after less than a dozen iterations, the program converges 
to a stable disk configuration. The problem we are solving 
is intrinsically two dimensional, but our MC code is based 
on cubic cells that are three dimensional, therefore we label 
such computations 2D disk models hereafter. 

Particular attention has to be paid to a good set-up 
of the grid used in the MC code. In the positive (x, y, z) 
directions, we use a basic grid of (350,350,50) cubes. At 
locations where the radiation field varies strongly, cubes 
are sampled into subcubes. Large gradients of the radiation 
field are expected close to the surface, within the extinction 
layer or at the inner wall of the disk. Estimates of the height 
of the extinction layer zq and its thickness £ ~ H/2 is given 
by Siebenmorgen & Kriigel (2010). 




1 1- , , 

10 100 
wavelength (um) 




radius (AU) 



Fig. 1. Top: IR emission of a T Tau disk viewed at ~ 45°, 
which is either computed using the 1+1D model (label 'ID') 
in a flat (green full line), a flaring (dashed green line) disk 
configuration, or the 2D disk. Iterations 1 to 9 of Eq. Q 
of the 2D disk are shown. Bottom: Midplane temperatures 
of the 1+1D flat disk and the various iterations of the 2D 
disk. 



The extinction layer is sampled with cubes having ry J5 
0.6. This is achieved by sampling cubes located close to 
the extinction layer in up to 20 subcubes. In total the disk 
has 10 7 cells, which provides a sampling with high spatial 
resolution. In each MC run a number of 10 8 photon packets 
are emitted from the star. The fairly high number of packets 
provide a high-energy resolution and is required to arrive at 
a stable estimate of the temperatures T(r, z), in particular 
near the midplane. 

4.3. Comparison: 1+1D versus 2D 

The SED and the midplane temperatures of a flat and flar- 
ing 1+1D disk are compared with those of a 2D disk model. 
As an example, we take a T Tauri star with parameters as 
specified in Tabled] For ease of comparison, the distribution 
of the dust column density S(r) is chosen to be identical in 
all three models. We take the one computed by the 2D disk 
as evaporation radius. The SED is compared in the upper 
panel of Fig. [TJ In a 1+1D disk the heating is proportional 
to the grazing angle. Since flaring disks have a larger graz- 
ing angle than flat disks they show higher dust re-emission 
luminosities and appear warmer. The assumption of 1+1D 
disk models is that stellar photons impinge upon the disk 
at a given radius with constant grazing angle. Such a sim- 
ple picture fails close to the dust evaporation zone. In early 
1+1D models the inner disk region was treated as an op- 
tically thick vertical wall at constant temperature (Natta 
et al. 2001; Dullemond et al. 2001). Then the dependency 
of gas pressure on the dust evaporation temperature is in- 
cluded, resulting in a curved evaporation zone (Isella & 
Natta 2005; Kama et al. 2009). In the MC scheme the ra- 
diation transport in radial direction is included, which, be- 
cause of the fairly difficult geometry, is otherwise not easy 
to implement in ray-tracing techniques as used in 1+1D 
disks. The SED of a 2D disk is very distinct with a much 
stronger mid-IR emission and overall natter appearance. 
We also show in Fig. [T] the convergence of the SED and 
the midplane temperature of the MC model as computed 
during the first nine iterations of the density structure fol- 
lowing Eq|4l The midplane temperature distributions of the 
1+1D disks are described well by a power law, whereas in 
the 2D disk striking up and downs with radius are derived. 
We extensively tested that the structure is not due to sys- 
tematic effects that could be introduced by the noise caused 
by a pure photon statistics or inadequate sampling of the 
MC cells. Following Eq. [3] the midplane temperature is di- 
rectly linked to the density, and one expects a hilly path, 
as derived in the lower panel of Fig. [I] to be reflected in 
the overall observed disk structure (Sect. [5]). 

4.4. SED dependency on dust opacities 

A realistic description of dust in protoplanetary disks is an 
open and controversial issue, as are the applied dust ab- 
sorption and scattering efficiencies. Because these parame- 
ters directly enter into the radiative transfer equation and 
the MC method, we wish to study their impact on the SED. 
As an example we take a T Tauri star with the parameters 
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Fig. 2. IR emission of 2D disks viewed nearly face-on 
20°) computed by applying different dust opacities. 



in Table Q] Three different dust models are distinguished: 
a) ISM like grains, b) fluffy grains with particle radii up to 
a + = 0.3/im, and c) fluffy grains with a + = 33/im. Other 
parameters remain as given in Sect. [2] The SEDs of the 
disks are compared using identical distributions of the col- 
umn density £(r). Since the dust models have different vi- 
sual extinction coefficients Ky, the vertical optical depth 
at 1AU of 10 4 for our standard dust (model c) needs to be 
increased by a factor 7.7 for ISM dust and by a factor 8.5 
for fluffy grains with size distribution, as for ISM grains. 
SEDs of the 2D disks are displayed in Fig. [3J where one 
notices that: 

i) The flux in the far-IR peaks at longer wavelength 
when applying fluffy grains. 

ii) The silicate emission profile depends on grain poros- 
ity, particles size, and dust temperature (Voshchinikov & 
Henning, 2008). In the band the cross section peaks at 
9.5/im for ISM dust and at 9.8pm for the fluffy composites. 
The fcature-to-continuum ratio of the extinction cross sec- 
tion in the 10/im band, when averaged over the dust size 
distribution, is 25% higher for ISM than for fluffy grains 
of same size. Nevertheless, because of the detailed temper- 
ature structure, model b), which is made of fluffy grains 
with sizes typical of the ISM, shows the most striking sili- 
cate emission feature. 

iii) The disk with fluffy grains up to a + = 0.3/im have 
a similar submillimeter slope as the one with 100 times 
larger, a+ = 33/im, particles (Fig. 0). This agrees with 
observations of protoplanetary disks where the silicate fea- 
ture is not correlated with the millimeter slope of the SEDs 
(Lommen et al., 2010). 

4.5. Disks plus halos 

Protoplanetary disks are thought to be formed from their 
initial spheroidal configuration of the protostellar envelopes 
(Watson et al., 2007). It is speculated that there is some 
residual dust left or that it is replenished at high altitudes 
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Fig. 3. IR emission of 2D disks viewed nearly edge-on (~ 
20°) with (dashed) and without (full line) a dusty halo. The 
disks arc cither heated by a T Tauri star (top) or a Herbig 
Ae star (bottom). 



above the disk (Vinkovic et al., 2006). Dust in a halo de- 
cays through radiation pressure and gravity and may be 
replenished by dynamical processes, such as outflows and 
winds from the disk surface, or by planets in highly ec- 
centric orbits (Krijt & Dominik, 2011). Miroshnichenko et 
al. (1999) added the emission of an optical thin halo onto 
the SEDs computed for optically thick disks. They found 
that the halo dominates the IR emission and provides ad- 
ditional disk heating on top of the direct stellar radiation. 
The MC scheme allows a self-consistent treatment of the 
radiative transfer in a configuration of a disk plus halo. In 
the T Tauri disk models (Table Q} we add a halo so that the 
minimum dust density of Eqs. [3] and 0] in each cell of the 
MC grid exceeds p(r, z) > p ha i = 1.5 x 10~ 18 (g-dust/cm 3 ). 
This gives an optical depth for the halo of ry ~ 0.4. The 
resulting SED of a disk plus halo is shown in Fig. [3] The 
SED appears warmer than the pure disk with phaio = 0. 
In the example, the disk plus halo model produces a fac- 
tor ~ 3 stronger IR peak emission and less near-IR flux 
for wavelengths below 6/im. The halo provides additional 
heating on top of the direct stellar light. Their midplanc 
temperature becomes a factor ~ 2 warmer than derived in 
a disk without halo within 1-10 AU of the star. 
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We repeat the exercise for the more luminous Herbig Ae 
stars with the parameters given in Table Q] With the same 
minimum dust density, /Ohaio, the optical depth of the halo 
is ry ~ 0.6 as for the T Tauri star. The resulting SEDs 
are shown in Fig. [3] The halo dominates the IR emission 
for wavelengths >3/im and produces warmer dust in the 
midplanc than in pure disks. The midplane temperature of 
the Herbig Ae stars are displayed in Fig. [5] The midplane 
temperature in the pure disk shows a rapid decline in the 
inner 3-4 AU, followed by a 1/r 6 dependency up to 12 AU 
and a strong decline farther out. The disk plus halo model 
is smoother and fit by T m id ~ 1/r in the 2-15 AU region 
(Fig.®. 



5. Ring structure of protoplanetary disks 

We notice that the midplane temperature in 2D disks of 
T Tauri stars show a wavy structure in particular near the 
zone where terrestrial planets are supposed to form (Fig. [3]). 
Here we investigate how this temperature behavior mani- 
fests itself in the appearance of the disk. First we study the 
height of the photosphere which we define as the height of 
the bottom of the extinction layer, Zbot, where the vertical 
optical depth tl(V) = 1. In FigH]we show the height of 
the disk photosphere as a function of radius, expressed as a 
unitless ratio Zbot/r. We do this for the T Tauri disk models 
with and without a halo using parameters of Table [T] For 
1+1D disks the midplane temperature is a smooth func- 
tion declining with radius, which can be approximated by 
a power law, and so is the scale height which smoothly in- 
creases with distance (Figj4]). The geometrical thickness of 
the 1+1D disk increases with radius. 

The 2D disk near the evaporation zone is puffed up, and 
then its thickness declines similar to the decrease in the 
midplane temperature. The temperature reduction is to be 
understand by shadowed region of the disk surface where 
light from the star is extincted and where grain heating 
becomes less efficient. The shadow is located close behind 
from the puffed-up inner rim. The disk becomes thicker 
with increasing distance, because the gravity decreases, so 
that there is a point where the shadow becomes less effi- 
cient, and the disk is exposed to direct stellar light. There 
again the disk is puffed up followed by a second shadow, 
which is located at about 1 AU in our example and a third 
one near 4AU. Farther out (>10 AU) the midplane tem- 
perature drops to below 30 K because the stellar heating of 
the passive disk is inefficient, and another distortion in the 
disk surface is not visible in most cases. The detailed struc- 
ture of passively heated T Tauri disks are shown in Fig0] 
within the region of terrestrial planets. If one considers ha- 
los, then gaps and ring-like structures are dimmed and the 
disk appears thicker because their midplane temperature 
is warmer. The disk surface is visualized in a three di- 
mensional representation by rotating the function Zbot (r) jr 
along the midplane (Fig. 2]). The innermost puffed-up rim, 
several gaps, and an overall wavy surface structure of the 
disk are visible. In Fig. [5] we display images and surface 
profiles for the mid-IR. At other wavelengths the emission 
structure is similar. In the image of the disk without halo 
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Fig. 4. Top: The relative height of the disk photosphere, 
Zbot/r, as a function of radius of a disk heated by a T Tauri 
star without (full line) and with (dashed) halo. The scale 
height over radius, H/r, of a flat 1+1D disk (dashed) is 
shown for comparison. Gaps between ring structures are in- 
dicated. A shaded representation of Zbot/f along the (x,y)- 
plane is displayed for a disk plus halo configuration (middle) 
and a pure disk (bottom) . Color bars give the relative gray 
scale (%) in units of Zb ot /r. 



there is a wide gap opening farther out and smaller gaps 
visible closer to the star. This structure is dimmed when a 
halo is considered. 
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Fig. 5. Mid-IR images of disks heated by a T Tauri star at 
distance of 50 pc and viewed at 30°. The image of a pure 
disk (top) and a disk plus halo is shown (middle). Color 
bars are in log(Jy/arcsec 2 ). Bottom: Corresponding surface 
brightness distribution measured as cut along S through 
origin of the image. Gaps between emission rings are indi- 
cated. 
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Fig. 6. Top: The midplane temperature and relative height 
of the disk photosphere, ZbotA, as a function of radius of a 
2D disk heated by a Herbig Ae star. Models are shown in a 
configuration of a pure disk (full line) and a disk plus halo 
(dashed). Power law fits (dotted) to the midplane temper- 
ature distribution are as labeled. Gaps between ring struc- 
tures are indicated. Bottom: A shaded representation of 
Zbot/r along the (x,y)-plane is shown for a pure disk model. 
Color bar gives relative gray scale (%) in units of ZbotA*- 
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Fig. 7. Top: Mid-IR image of a disk heated by a Herbig 
Ae star at distance of 50 pc and viewed at 30° . Color bar 
is in log(Jy/arcsec 2 ). Bottom: Mid-IR surface brightness 
distribution measured as a cut along S through the origin 
of an emission image of a pure disk (full line) and that 
computed in a disk plus halo configuration (dashed). 



Disks of the more luminous Herbig Ae stars have an 
evaporation zone farther out than around T Tauri stars. 
At these distances the gravitation is reduced and the disks 
become thicker than for T Tauri disks. The height of the 
Herbig Ae disks and a shadowed presentation is displayed 
in Fig. El When compared to T Tauri stars (Fig. H]) , the ra- 
tio Zbot(r)/r is indeed higher and has less structure. There 
is a pronounced puffed-up inner rim causing a far reaching 
shadow, so that a second one farther away than 30AU is 
only envisaged. However, this far away the midplane tem- 
perature has dropped to below 30 K (Fig. and other ef- 
fects than stellar heating may dominate the disk structure. 
We display the mid-IR image and surface profile of the disk 
heated by a Herbig Ae star in Fig. [7J At least one striking 
gap and a ring in the disk are identified. 



The infrared appearance of passive disks around low and 
medium mass stars has been discussed. Dust particles are 
made either of bare material or fluffy composites of silicate 
and carbon. Grain sizes range from 160A to 3000A as de- 
rived for the ISM and up to 33/im assuming grain growth 
in the disks. The disks are heated by stellar light with lumi- 
nosities between 2 L Q and 50 Lq . The disks have an inner 
radius that is roughly determined by the dust evaporation 
temperature of the grain material. The surface density of 
the disks was set close to what is estimated for the mini- 
mum mass of the solar nebula. Disks are configured with 
or without halos. Two hydrostatic and geometrically thin 
disk scenarios were examined for which we assume that gas 
and dust are mixed and have the same temperature. 

i) In 1+1D the disks were divided in small ring segments 
in which the radiation transfer was solved. For each ring, 
a slab geometry was applied, and it was assumed that in 
deeper layers the disk is isothermal. The transport of radi- 
ation in the radial direction was ignored. From the surface 
of 1+1D disks, the star is always visible. A flat or flaring 
structure of the disk was considered. This type of disk is in 
widespread use in the literature. 

ii) In 2D disks the hydrostatic and radiative transfer 
equations were solved in an iterative scheme. For the latter 
an accurate MC method was presented that can be ap- 
plied to arbitrary dust geometries. We used an adaptive 
three-dimensional Cartesian grid where cubes are divided 
into subcubes whenever required, for example, close to the 
disk surface. The MC method is vectorized and the paral- 
lelization is realized to work on graphics cards or conven- 
tional processing units and provides a speed-up roughly 
proportional to the number of processor cores or parallel 
threads, available. On our conventional computer, the vec- 
torized code is a factor 100 faster than the scalar version. In 
regions of very high optical depth, for example close to the 
midplane of the disk, photon packets may get trapped. In 
this case the algorithm is further accelerated by a modified 
random walk procedure. The 2D disks provide a more real- 
istic description than the 1+1D disks: In the RT solution of 
the 1+1 D disks it is assumed that the radial transport of ra- 
diation can be ignored and that layers with vertical optical 
of Ttop ^ 20 are isothermal. These assumptions break for ex- 
ample near the inner regions of the disk. On the other hand, 
in the 2D disks the radial transport of radiation is included 
and no assumptions are made that regions of the disks need 
to be isothermal. Therefore a self consistent treatment of 
the hydrostatic balance is only provided in the 2D disks. 
Our main findings are that: 



— 1+1D disks have a smooth surface structure, and their 
midplane temperature is approximated by a power law 
well. They gravely underestimate the mid IR emission 
when compared to 2D disks. 

— Halos substantially increase the IR emission and the 
temperature of the midplane, which is smoother than 
obtained in disks without halos. 
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— In 2D disks the midplane temperature shows up and 
downs with radius (FigJlJ. Such a hilly structure is also 
revealed in their surface structure. 

— Emission images of 2D disks display gaps and ring-like 
structures in particular in the region of terrestrial plan- 
ets. Disk structures are caused by pufhng up the disk 
surface and followed up by their shadows behind. The 
disk structure is dimmed when halos are considered. 
Rings and gaps are more pronounced in T Tauri stars 
than in Herbig Ae stars. 

The detected gap and ring structures of the disks are 
only caused by the assumption of hydrostatic equilibrium 
and the derived detailed vertical temperature profiles of the 
disk. On the other hand, it is easy to imagine, without the 
need of such detailed computations, that behind a puffed- 
up disk the heating is dimmed and therefore that the dust 
temperatures are lower. This causes that in those shadowed 
regions the height of the disk is smaller unless it is again 
exposed to direct stellar light. Therefore the detected gap 
and ring-like structures are a plausible effect of the shad- 
ows, but they have not yet been reported by others. We 
feel that this is important because ring-like structures of 
protoplanetary disks are often interpreted as "fingerprints" 
of the planet formation process, whereas in the models pre- 
sented planets have not been considered. 

Acknowledgements. We are grateful to Endrik Kriigel for discussions 
and helpful suggestions. 
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